#Created by Zhaohui Yang
#Modified by Jinchi Lu (jinlu@ucsd.edu)
#elastic pressure independent material
# plane strain,  single element,
# push over analysis (Applied load linearly increases with time),
# English units (in, s, lbf)
#
#         4     3
#         ------- --> F  (loads applied to node 3)
#         |     |
#         |     |
#         |     |
#        1-------2   (nodes 1 and 2 fixed)
#         ^     ^
wipe
# define UNITS ---------------------------------------------
# Modified from Silvia Mazzoni & Frank McKenna, 2006
set in 1.;              # define basic units -- output units
set lbf 1.;             # define basic units -- output units
set sec 1.;             # define basic units -- output units
set ft [expr 12.*$in];          # define engineering units
set ft2 [expr $ft*$ft];
set ft4 [expr $ft2*$ft2];
set kip [expr 1000.*$lbf];
set ksi [expr $kip/pow($in,2)];
set psi [expr $ksi/1000.];
set pcf [expr $lbf/pow($ft,3)];         # pounds per cubic foot
set psf [expr $lbf/pow($ft,2)];         # pounds per square foot
set in2 [expr $in*$in];                                 # inch^2
set in4 [expr $in*$in*$in*$in];                 # inch^4
set g [expr 32.2*$ft/pow($sec,2)];      # gravitational acceleration

set psfTokPa   4.788000e-002
set psiTokPa   6.894757e+000
set ksiTokPa   6.894757e+003
set ftTom   3.048000e-001
set inTom   2.540000e-002
set kipsTokN   4.448222e+000
set kipsftTokNm   1.355818e+000
set pcfToton_m3   1.601800e-002

set kPaTopsf  [expr 1./$psfTokPa];
set kPaTopsi  [expr 1./$psiTokPa];
set kPaToksi  [expr 1./$ksiTokPa];
set mToft [expr 1./$ftTom];
set mToin [expr 1./$inTom];
set kNTokips [expr 1./$kipsTokN];
set kNmTokipsft  [expr 1./$kipsftTokNm];
set ton_m3Topcf  [expr 1./$pcfToton_m3];
set kkg_m3Toslug_ft3 1.94055


#
#some user defined variables
#
set massDen  [expr 2.*$kkg_m3Toslug_ft3/$ft4];       ;# solid mass density
set fluidDen [expr 1.*$kkg_m3Toslug_ft3/$ft4];       ;# fluid mass density
set cohesion [expr 30*$kPaToksi*$ksi]         ;
set peakShearStrain 0.1 ;
set E1      [expr 90000.0*$kPaToksi*$ksi]      ;#Young's modulus
set poisson1 0.40 ;
set G [expr $E1/(2*(1+$poisson1))] ;
set B [expr $E1/(3*(1-2*$poisson1))] ;
set press   [expr 0.*$kPaToksi*$ksi]           ;# isotropic consolidation pressure on quad element(s)
set period   2          ;# Period of applied sinusoidal load
set deltaT   0.01       ;# time step for analysis
set numSteps 2000       ;# Number of analysis steps
set gamma    0.5        ;# Newmark integration parameter
set pi 3.1415926535     ;
set inclination 0       ;
set unitWeightX [expr  ($massDen-$fluidDen)*$g*sin($inclination/180.0*$pi)] ;# buoyant unit weight in X direction
set unitWeightY [expr -($massDen-$fluidDen)*$g*cos($inclination/180.0*$pi)] ;# buoyant unit weight in Y direction
set loadIncr [expr 1.*$kNTokips*$kip]          ;# Static shear load

set frictionAng 0
set refPress  [expr 100*$kPaToksi*$ksi];
set pressDependCoe 0.
set numberOfYieldSurf 20


####################################################################
#create the ModelBuilder
model basic -ndm 2 -ndf 2

# define material and properties
nDMaterial PressureIndependMultiYield 2 2 $massDen $G $B  $cohesion $peakShearStrain $frictionAng $refPress $pressDependCoe $numberOfYieldSurf
#nDMaterial FluidSolidPorous 1 2 2 2.2e6

# define the nodes
node 1   0.0 0.0
node 2   [expr 1.0*$mToft*$ft] 0.0
node 3   [expr 1.0*$mToft*$ft] [expr 1.0*$mToft*$ft]
node 4   0.0 [expr 1.0*$mToft*$ft]

# define the element      thick  material      maTag   press  density  gravity
element quad  1  1 2 3 4  [expr 1.0*$mToft*$ft]   "PlaneStrain"     2   $press  0.0    $unitWeightX  $unitWeightY

updateMaterialStage -material 2 -stage 0

# fix the base in vertical direction
fix 1 1 1
fix 2 1 1

#tie nodes 3 and 4
equalDOF 3 4   1 2

#############################################################
# GRAVITY APPLICATION (elastic behavior)
# create the SOE, ConstraintHandler, Integrator, Algorithm and Numberer
system ProfileSPD
test NormDispIncr 1.e-12 25 0
constraints Transformation
integrator LoadControl 1 1 1 1
algorithm Newton
numberer RCM
analysis Static
analyze 2



#############################################################
# NOW APPLY LOADING SEQUENCE AND ANALYZE

# rezero time
setTime 0.0
wipeAnalysis

# create a LoadPattern with a Linear time series
pattern Plain 1 Linear {
    load 3 $loadIncr 0.0    ;#load applied in x direction
}

# create the Analysis
constraints Transformation ;  # Penalty 1.0e18 1.0e18  ;#
test NormDispIncr 1.e-12 25 0
algorithm Newton
numberer RCM
system ProfileSPD
rayleigh 0.0 0.0 0.000 0.
integrator Newmark $gamma  [expr pow($gamma+0.5, 2)/4]
analysis VariableTransient



#create the recorder
recorder Node -file disp.out   -time  -node 1 2 3 4 -dof 1 2 -dT  0.01 disp
recorder Node -file acce.out  -time  -node 1 2 3 4 -dof 1 2 -dT 0.01 accel
recorder Element -ele 1 -time -file stress1.out -dT 0.01 material 1 stress
recorder Element -ele 1 -time -file strain1.out -dT 0.01 material 1 strain
recorder Element -ele 1 -time -file stress3.out -dT 0.01 material 3 stress
recorder Element -ele 1 -time -file strain3.out -dT 0.01 material 3 strain

#analyze
set startT [clock seconds]
analyze $numSteps $deltaT [expr $deltaT/100] $deltaT 10
set endT [clock seconds]
puts "Execution time: [expr $endT-$startT] seconds."

wipe  #flush ouput stream